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Abstract. - We derive and solve a differential equation satisfied by the probability distribution 
of the work done on a single biomolecule in a mechanical unzipping experiment. The unzipping 
is described as a thermally activated escape process in an energy landscape. The Jarzynski 
equality is recovered as an identity, independent of the pulling protocol. This approach allows 
one to evaluate easily, by numerical integration, the work distribution, once a few parameters 
of the energy landscape are known. 



Introduction. - The introduction of micromanipulation techniques has dramatically im- 
proved our knowledge of physical and chemical properties of biological molecules. Such tech- 
niques have been used to probe the structure of proteins [1-5] and nucleic acids [6] . A typical 
experiment consists of pulling the free end of a biomolecule with a controlled force, while its 
end-to-end distance is measured at the same time. It has been suggested that the study of 
the kinetics of bond breakage under different loading rates can provide many informations on 
the molecule internal structure, and in particular allows one to measure the strength of the 
molecular bonds, and to associate to them a position along the molecular structure [7]. The 
loading-rate dependent kinetics experiments on biomolecules have been interpreted in terms 
of thermally activated escape from bond states over a succession of energy barriers, along a 
one-dimensional energy landscape [8-11]. 

Usually, because of technical limitations, the molecule pulling process is characterized 
by time scales much faster than the typical molecular relaxation time. This prevents the 
possibility to perform the experiment in quasiequilibrium conditions and thus to obtain direct 
measurements of the thermodynamic state variables. This difficulty can be overcome by 
exploiting the remarkable equality derived by Jarzynski [12,13] (Jarzynski equality, JE), and 
extended by Crooks [14] , which allows one to obtain the free energy difference AF between two 
equilibrium states by evaluating the average of exp(— f3W), where W is the work performed 
on the system during the thermodynamic transformation: 

e -0AF = ^-W} . (1) 
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Here (3 = and the average on the rhs is performed over all possible realizations of the 

process starting from the equilibrium ensemble at a given temperature T. Since this equality 
holds in general, the information of AF can be gathered also by processes so fast that they do 
not leave the system at equilibrium. On the other hand, the process must be sampled a large 
number of times in order to obtain a reliable estimate of AF, as discussed in [15, 16]. In the 
same references it is argued that the number of pulling experiments needed to achieve good 
statistics in the estimate of AF increases as the transformation are made more irreversible, 
for example as the pulling rate is increased. Furthermore the JE gives no information on the 
probability distribution of work: while the JE magnifies the rare trajectories with W < AF, 
a direct measurements of the probability distribution for these value of W is rather difficult. 
On the other hand, the effect of such trajectories on the average thermodynamic variables 
might become significant when the energies involved range between a few and tens of fc^T's, 
as in the case of molecular bonds. 

The aim of this paper is thus to provide an effective method to evaluate the work probability 
distribution of a molecular pulling process, which can be described as an escape process in an 
energy landscape, once the main features of the landscape are known. Such a method can be 
applied independently of the irreversibility of the process, i.e., of the pulling protocol, and for 
any value of the maximum force. 

The article is organized as follows. We first describe our model of escape process. We then 
introduce a set of differential equations describing the time evolution of the work distribution 
probability <p(W, t) and show that the JE follows as an identity. We next consider a simple 
case of escape process, and discuss the behaviour of the work distribution probability 4>(W, t) 
obtained for it. 

The model. - In a typical unzipping experiment, a force f(t) is applied on one end of the 
biomolecule, and its elongation x(t) is monitored. Within some limits, x(t) can be considered 
as a collective coordinate for the system, spanning a one-dimensional (free) energy landscape 
E(x). This energy landscape will in general be characterized by a set of N minima of energy e% 
at position Xi, with i = 0, . . . , N — 1, and by a set of N— 1 maxima of energy Ej at position Xj 
with j = 1, . . . , N — 1, see figure ^ The escape from one minimum e, over the next maximum 
Ei+i can be viewed as the breaking of a given molecular bond. The evolution of the collective 




Figure 1 - Generic energy landscape, representing the succession of bonds in a biomolecule. 

coordinate will be governed by a stochastic process satisfying detailed balance. If however 
the variations of the energy E between maxima and minima are large enough with respect 
to kftT, it is possible to describe this stochastic process as a Markov process with discrete 
states, corresponding to the extremals of the energy, while the transition probabilities will 
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be given by the Kramers expression. For simplicity, we shall now consider only transitions 
between nearest neighbour energy minima j — > i, i.e. i = j ± 1. However, we will argue that 
the results we obtain in the present paper are rather general and can be extended to the case 
of transitions between any pair of states j and i. Thus the transition rate from one minimum 
j to a neighboring one i, over the corresponding energy barrier, is given by 

hj(t) = ljq exp {-(3 [Eij - ej - f(t)(X i:j - Xj)}} , (2) 

where the energy Eij and the position of the barriers are respectively given by 

Eij £'max{!j} 5 Xij ^max{i,j} ) (3) 

and where luq is some attempt rate whose value depends on the system characteristics. Let 
Pi(t) be the probability that the system is in state i at time t: using a Kramers formalism the 
time evolution of these probabilities is described by the following set of differential equations 

= E M*)pj(*) - **(*)*(*)] ■ (4) 

j 

The jump process, described by the set of equations is clearly Markovian and it preserves 
the starting equilibrium ensemble for any fixed value of /. These conditions, together with 
the assumption that the energies involved in the process are finite, are sufficient for the JE 
to be recovered for the dynamics described by eqs. (0J, as shown by Jarzynski [12, 13] and 
Crooks [14]. 

In the following Z t will denote the the partition function of the canonical ensemble at 
temperature T and with an applied force equal to fit): 



E 



exp{-/?[e;-/(i)x;]} . (5) 



Work probability distribution. - The probabilities pi (t) introduced in the previous section, 
are not sufficient to describe the probability distribution of the work done on the system. We 
are thus interested in the joint probability distribution 4n(W, t) that the system is in state i 
at time t, while the total work done on it is equal to W . 

If the system is in the state i at time t, the work done SW in the time interval St is given 

by H 

SW = -5t^L Xt . (6) 

By expanding to first order in St and SW, the probability distribution <fii(W + SW,t + St), 
it is easy to verify that this function satisfies the set of differential Chapman-Kolmogorov 
equations: 

- E Mttow*) - m«)*w*)] + h*)«y$m . (7 ) 



( 1 )In the present paper we use the ensemble where the force is the externally controlled parameter, while the 
extension of the molecule Xi fluctuates. This implies that the effective work done on the system, as the force 
increases of Sf, is SW = —XiSf. In pulling experiments where an optical tweezer is used to pull one of the free 
ends of the molecule, the actual controlled parameter is the position of the focus of the optical trap. In this 
situation both the applied force and the molecule extension fluctuate. Here we assume that the fluctuations 
in fit) are small, even if the conditions under which this assumption is valid deserve further investigation. 
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The probability Pi(t) and the probability distribution cfii(W,t) are connected by the relation 

dWMW,t)=Pi(t), (8) 



provided that 4>i(W, t) satisfies the additional boundary conditions <^(±oo, t) — 0. By in- 
tegrating both sides of equation J7J with respect to W, one recovers the set of differential 
equations Q. 

We want now to show that the JE is satisfied for a system characterized by a generic 
energy landscape as that represented in fig. ^ whose work probability distributions evolve 
according to eq. J7J). Let <f> l j(W,t) denote the probability distribution of the work in the state 
j at time t, when the system initial state is I. They satisfy the set of equations 10 with the 
initial conditions 

S(W), if j = i; 

0, if j^l; 



d(w,t = o) = 



(9) 



Thus the average of e & w up to time r is given by 



(e^l = E 



exp(-/?e ; ) 



-pw 



dW = 



where 



- pW (j){W,T)dW , (10) 



(11) 



is the total work distribution probability. Our goal is to show that the rhs of eq. p0(l is equal 
to exp(— (3AF). Let us define the quantity Mi(t) as 



M i {t)= e-^Y. 
J i 

the rhs of eq. (|10f) can thus be written as 



-0e t 



-4(w,t)dw , 



By taking the time derivative of Mi(t), and substituting eqs. 0) and (|12fl . we obtain 
dMi(t) 



0W Y, e ^lti{w,t)dw 



dt J ^ Z dt 



(12) 



(13) 



(14) 



= ]T[MWW-Mi)W)] 



-pw 



-pe, 



E -^-f(t>iM\(W,t)dW .(15) 



The last term in eq. (jT3|l can be integrated by parts, and using the boundary condition that 
4>\(W,t) goes to zero faster than exp((3W) as W — > -co (which follows from the existence of 
the average in eq. Q), we find that Mi(t) satisfies 



= E NW m jW - MWt)] + Pf{t) Xt Mr{t) , 



(16) 
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with the initial conditions Mj(0) = e~^ ei /Zq, which follow from the initial value conditions 
for the functions <j>j(W,t), eq. 10, and from the definition of M;(t), eq. I|12[). It is easy to 
verify that the functions 

Mi(t) = «P[-^-/(*)^)l (17) 
Zq 

are the solutions of the set of equations (|16|l with the corresponding initial conditions. Sub- 
stituting this last result into eq. I|13|l . we finally find 

/ e-w m T )dW = Y «P[-^-/(^)l = ^ y (18) 
J Z Z 

which verifies the JE. 

Work distribution probability: a case study. We consider now a simple system with 
three minima and two maxima in the energy landscape. The values of the energy minima 
and maxima, expressed in fcgT units, and their positions, expressed in nm, are {eo = 0, xq = 
0}, {E 1 = 10, X 1 = 0.6}, {ei = 6, xi = 0.8}, {E 2 = 16, X 2 = 1.8}, {e 2 = 12, x 2 = 2}. 

We consider here the case of pulling processes where the force increases linearly with time, 
i.e. / = r ■ t, where r is the pulling rate. In order to simplify the notation, let 



e 



-pe, 



^t) = —^(W,t) . (19) 
The work probability distribution (f)(W, t) is thus given by 

c)>{W,t) =Y,${W,t) . (20) 

We obtain the probability <p(W,t), for two values of the pulling rate r = 1, 10 pN/s, by 
solving the set of equations J7J) with the parameter choice as indicated above, taking the 
initial value conditions as given by eq. 10 and the attempt frequency u>q = 4.4 s _1 . With this 
choice of parameters, we obtain the zero-force transition rate fco = i^ocxp [— /3(Ei — ej-i)] = 
2 x 10 -4 s _1 , in agreement with the zero-force transition rate found in the mechanical unfolding 
of a simple RNA molecule [6] . 

The results at different times arc shown in figure [3 for r = 1 pN/s and r = 10 pN/s. In 
order to analyze the behavior of the tails of the work distribution, the quantity <fi(W, t) for 
r = 1 pN/s is plotted in a log- linear form in figure |3 for the smallest and the largest times 
considered. Recalling that <f>(W,t) is the sum of the work probability distributions <j> l j(W, t) 
along all the possible trajectories, see eq. I|2U|). in the same figures the main contributions to 
the total distribution probability of the work are plotted. Inspection of figure |21 suggests that 
at small times and for small values of \W\ the work probability distribution is dominated by 
those trajectories which start and finish in the state xq. On the other hand, at any time, at 
large \W\, 4>(W, t) is dominated by those trajectories which finish in the rightmost state x 2 , 
as expected. In particular, in the very large \W\ regime, the distribution <fi 2 (W,t) determines 
the behavior of <fi(W, t). 

Thus, at any time, a given distribution function <fij(W,t) dominates the behavior of the 
total work distribution probability cf>(W, t) in a given range of W. In each of such work ranges, 
the function (f)(W, t) can be well fitted by a gaussian distribution (fits not shown), i.e., the 
function (f>(W, t) is a superposition of several gaussian functions, each with different mean 
and variance, rather than a single gaussian function. The gaussian distribution of the work, 
expected for a slowly perturbed system as discussed in [16, 17], is thus not recovered here for 
the pulling rate r = 1 pN/s, which is a lower limit in the pulling experiments of biomolecules. 
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Figure 2 - Plot of the work distribution probability 4>(W, t), as defined by eq. 12i)l as a function of 
the work W at different times. Left: r — 1 pN/s. Right: r = 10 pN/s. Inset: 4>(W,t) at the smallest 
time here shown (left: t — 3.3 s; right: t = .33 s). 



Discussion. - In the present paper we have introduced a theory which allows us to obtain 
the work probability distribution performed on a biomolccules during a pulling experiment. 
The pulling process has been described as a jump process over a succession of energy barriers. 
This approach allow us to describe the time evolution of the population of the minimal energy 
states using a set of stochastic differential equations. By discretizing such process, we show 
that the probability of the discrete trajectories satisfies the Jarzynski equality independently 
of the pulling protocol. The derivation of the JE within the current formalism can be viewed 
as the discrete-state-space analogue of the analysis in Sec. II of ref. [13] and in Sec. IV of 
ref. [18]. 

We then derive the set of differential equations describing the time evolution of the work 
probability distributions for each trajectory connecting two arbitrary initial and final states of 
the system. By summing up these distribution functions over all the initial and final states, we 
obtain the total work probability distribution cf>(W,t). We show that the JE is recovered for 
such a distribution function: in this sense, the average over an infinite number of trajectories 
of eq. (JIJ, is replaced by an average over the probability distribution cf>(W,t) as in eq. (|ll)fl . 
In deriving this results, we have assumed that the main features of the energy landscape, i.e., 
the height and the positions of the wells and of the barriers, are known. For real biomolccules 
these quantities have to be measured independently for the function 4>{W, t) to be numerically 
calculated. By considering a very simple energy landscape, we calculate numerically the 
work probability distribution, and find that it exhibits a non-gaussian behaviour, being rather 
a superposition of gaussian functions, each corresponding to distinct trajectories between 
different initial and final state. 

We point out that our results still hold if we consider a jump process where transitions 
between any pair of states are allowed. In fact, as long as the transition rates kij satisfy the 
detailed balance condition, eqs. (|14I18|I ar e still verified, and thus one recovers the results 
obtained for the case of jumps between successive states. 



* * * 



We are grateful to F. Ritort for introducing us to the topic and for interesting discussions. 
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Figure 3 - Log-linear plot of the work distribution probability <f>(W,t) as a function of W for r = 
1 pN/s, and of the main contributions to the sum on the right hand side of eq. (I20L at short times 
(ti = 3.3) and at long times (£2 = 63 s). 
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